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Abstract 



The electron transfer kinetics of mixed- valence systems is studied via solv- 
ing the eigen-structure of the two-state non-adiabatic diffusion operator for a 
wide range of electronic coupling constants and energy bias constants. The 
calculated spectral structure consists of three branches in the eigen-diagram, a 
real branch corresponding to exponential or multi-exponential decay and two 
symmetric branches corresponding to population oscillations between donor 
and acceptor states. The observed electronic coherence is shown as a result of 
underdamped Rabi oscillations in an overdamped solvent environment. The 
time-evolution of electron population is calculated by applying the propagator 
constructed from the eigen-solution to the non-equilibrium initial preparation, 
and it agrees perfectly with the result of a direct numerical propagation of the 
density matrix. The resulting population dynamics confirms that increasing 
the energy bias destroys electronic coherence. 
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I. INTRODUCTION 



Quantum coherence in the dynamics of condensed phase systems has become a subject of 
recent experimental and theoretical studies. A central issue is the observability of electronic 
coherence in electron transfer systems given the fast dephasing time in many-body quantum 
systems. Experimentally, with the advance in ultrafast laser technology, oscillations in elec- 
tronic dynamics have been observed in photo-synthetic reaction centers and other electron 
transfer systems and are believed to arise from vibrational and/or electronic coherence. 0~i 
Accurate measurements on photo-induced electron transfer in mixed-valence compounds 
have demonstrated oscillations in electronic populations on the femtosecond time-scale.Bi 
Theoretically, detailed path-integral simulations suggest that such oscillations take place in 
electron transfer systems with large electronic coupling constants and are sensitive to the 
initial preparation of the bath modes associated with the transfer processes. Lucke et a/.i 
extended the non-interacting blip approximation to incorporate the non-equilibrium initial 
preparation and carried out extensive path-integral quantum dynamics simulations for elec- 
tron transfer reactions. According to their findings, large-amplitude oscillations are most 
likely to be observed in symmetric mixed-valence systems that are nearly adiabatic and 
with initial configurations that are centered in the Landau-Zener crossing region. Using the 
transfer matrix technique]! Evans, Nitzan, and Ratneri calculated short-time evolution for 
the photo-induced electron transfer reaction in (NH3)5Fe n (CN)Ru m (CN) 5 . Their results 
show fast oscillations in the electronic population on the short time-scale (20 fs) followed 
by a slower population relaxation on the long time-scale (100 fs). They pointed out that 
these fast oscillations arise as the wave-function oscillates coherently between the donor 
and acceptor states. The calculated long-time decay rate is considerably smaller than the 
prediction by the golden-rule formulae,!!! confirming the inadequacy of non-adiabatic rate 
theory in studying mixed- valence systems. 



In fact, a simple classical argument helps understand the nature of the observed os- 
cillations. As a function of the ratio between A (the bath reorganization energy) and V 
(the electronic coupling constant), there is a thermodynamic transition from the localized 
electronic state in a double-well potential to the delocalized electronic state in a single well 
potential.0Hll (i) In the localized regime (A > F), the large reorganization energy destroys 
electronic coherence; hence, electron transfer is an incoherent rate process, which can be de- 
scribed by the non-interacting blip approximation or golden-rule rate in the non-adiabatic 
limit and by transition state theory in the adiabatic limit (ii) In the delocalized regime 
(A < V), the electronic wave function extends to both the donor and acceptor states and 
electronic coherence persists over several oscillations.0 For mixed-valence compounds, the 
electronic coupling constant is estimated to be in the range of 10 3 cm _1 , which is in the same 
order as the reorganization energy.0'1 Therefore, the observed oscillations and relaxation in 
mixed-valence systems are the consequence of a highly non-equilibrium coherence transfer 
process. 

Due to the derealization nature of electronic states, an adiabatic pictureEl is more useful 
than the diabatic representation for analyzing the short-time dynamics in strongly-coupled 
systems. In this picture, electronic coherence arises from Rabi oscillations between two adi- 
abatic surfaces and decays because of electronic dephasing. Further, initial preparation and 
wave-packet dynamics can modulate Rabi oscillations and the overall electronic dynamics. 
Thus, the adiabatic representation provides a simple picture for mixed-valence systems as 
well as a simple analytical method to model fast electron dynamics initiated by laser pulses. 

As a general approach to describe condensed phase dynamics, we recently proposed a 
spectral analysis methodjl which is based on eigen-structures of dissipative systems instead 
of dynamic trajectories. An important application of the approach is to analyze a set of 
two-state diffusion equations, which was first used by Zusman to treat solvent effects on 
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electron transfer in the non-adiabatic limit. The analysis allows us to characterize multiple 
time-scales in electron transfer processes including vibrational relaxation, electronic coher- 
ence, activated curve crossing or barrier crossing. With this unified approach, the observed 
rate behavior, bi-exponential and multi-exponential decay, and population oscillations are 
different components of the same kinetic spectrum. Thus, several existing theoretical mod- 
els, developed for limited cases of electron transfer, can be analyzed, tested, and extended. 
In particular, rate constants extracted from the analysis bridge smoothly between the adia- 
batic and non-adiabatic limits, and the kinetic spectrum in the large coupling regime reveals 
the nature of the localization-delocalization transition as the consequence of two competing 
mechanisms. 

In this paper, the spectral analysis approach developed in Ref . |1^ is employed to study the 
electron transfer dynamics in mixed-valence systems. We invoke the non-adiabatic diffusion 
equation proposed by Zusman to describe the electron transfer process in the over-damped 
solvent regime. As discussed earlier, electron transfer in mixed-valence systems takes place 
in a different kinetic regime from the thermal activated regime described by Marcus theory. 
Thus, the time-scale separation is not satisfied, and multi-exponential decay and oscillations 
are intrinsic nature of electron transfer kinetics. As a result, the kinetic spectra exhibit bifur- 
cation, coalescence, and other complicated patterns. Careful examination of these patterns 
reveals the underlying mechanisms in mixed- valence systems. 

The rest of the paper is organized as follows: The spectral structure of the non-adiabatic 
diffusion equation is formulated in Sec. II. Numerical examples of the spectral structure 
of strongly mixed electron transfer systems are presented and discussed in Sec. Ill and 
concluding remarks are given in Sec. IV. 
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II. THEORY 



There have been extensive studies of the solvent effect on electron transfer dynamics 
in literature with various approachesJlHEl One of the most extensively studied models for 
quantum dissipation is the spin-boson HamiltonianJll'H 
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(1) 



where e is the energy bias between the two electronic states, V is the electronic coupling 
constant, a z and a x are the usual Pauli matrices, and {x a ,p a } represents the bath degree of 
freedom with mass m a , frequency ui a , and the coupling constant c a . In this model effects 
of the bath modes on the dynamics of the system can be described via the spectral density 
defined by, 
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^) = oE-7%-^)' (2) 

Equivalently, the spin-boson Hamiltonian in Eq. ([!]) can be separated into the electronic 
two- level part Htls an d the nuclear bath part Hb, 

Hsb = H TLS + H B . (3) 

The two-level part of the Hamiltonian can be explicitly written as 

Htls(E) = + U 2 {E)\2)(2\ + V{\1)(2\ + |2)(1|), (4) 

where the diabatic energy surfaces U\(E) and U2(E) are functions of the stochastic variable 
E, which represents the polarization energy for a given solvent configuration.!!] The trans- 
formation from the spin-boson Hamiltonian to the two-level system Hamiltonian has been 
shown in the literatureil'il by the identity, 

E({x a }) = J2c a x a . (5) 



It is worthwhile to mention that the polarization energy E was recognized as the reaction 
coordinate by Marcus in formulating non-adiabatic electron transfer theory.El Since the 
electron transfer process involves the collective motion of a large number of solvent degrees 
of freedom and the two-level system is linearly coupled to the harmonic bath modes in 
the spin-boson Hamiltonian in Eq. ([!]), the functional form for the free energy surface is 
harmonic,!! thus giving 

U 2 (E)= { -^^ + e , (7) 
where A is the reorganization energy, which is related to the parameters in Eq. ([[]), 

Considering the fact that electron transfer processes are usually probed at room temper- 
ature in polar solvents, we can treat the bath degrees of freedom in Hb classically. Then, 
the spin-boson Hamiltonian in Eq. (^) can be used to derive a two-level classical equation 
of motion, 

i^p(t) = Cp(t) = (C B + C TLS )p(t), (9) 

where iCs = {Hb, } is the Poisson operator for the classical bath and Ctls = [Htls, ]/h 
is the Liouville operator for the two level system. Explicitly, we express Eq. (|9|) in terms of 
the density matrix elements, 

pt = Lxpx + iV(pi 2 - P21), (10a) 

Pi = C2P2 - iV(pn - P21), (10b) 

P12 = C\2p\2 - iunpv2 + iV(p! - p 2 ), (10c) 

P21 = C21P21 + iuviPn - iV(pi - p 2 ), (10d) 
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where the Planck constant Ti is set to unity for simplicity, pi is the diagonal matrix element 
for electronic population, and is the off-diagonal matrix element for electronic coherence. 
Here, £ describes the relaxation process of classical bath, with Li defined on the free energy 
surface for the ith electronic state, and with £i 2 and £21 defined on the averaged free 
energy surface. This set of semi-classical two-state equations has been previously derived in 
different context by several authorsBlS It should be mentioned that the mapping from the 
spin-boson Hamiltonian into the Zusman model requires the Lorentzian form of the spectral 
density, 

J{u) = 2\^- 2 . (11) 

Furthermore, we note that many chemically and biologically important electron transfer 
processes take place in the over-damped solvent environment. Therefore, to describe the 
density matrix evolution in the electron transfer kinetics in the mixed-valence system, we 
invoke the non-adiabatic diffusion equation proposed by Zusman.E] Then, the bath relaxation 
operators in Eq. (||) are one-dimensional Fokker-Planck operators Cij, 

L-D 9 ( 9 +8 m{E) ) (12) 
Ci ~ DE dE [dE + P-d^ ) > (12) 

Cl2 ~ Al " 2 " ° E dE [dE +(3 ^E- ■ (13) 



where (3 = l/ksT, U and U12 are the average and the difference of the two free energy 
surfaces, respectively, 

m = ftw+ftw , (14) 

u 12 {E) = U l {E)-U 2 {E). (15) 
The energy diffusion constant De is defined as 

D E = ^A|, (16) 



where A% is the mean square fluctuation of the solvent polarization energy 

A| = (E 2 ) = 2Xk B T, 

and td = is the the characteristic timescale of the Debye solvent. The correlation 
function of the solvent polarization energy is given by 

C(t) = (E(t)E(0)) = A| exp(-nt). (17) 

Note that since the nuclear dynamics is modeled by the Fokker-Planck operator, the possi- 
bility of the vibrational coherence is excluded in this model of electron transfer dynamics. It 
is worthwhile to mention that one can obtain the non-adiabatic diffusion equation starting 
from the spin-boson Hamiltonian, by first deriving the evolution equation for the quantum 
dissipative dynamics, and then taking the semi-classical limit using the Wigner distribution 
functions, and finally assuming the over-damped diffusion limit .0 

We investigate the spectral structure of the non-adiabatic diffusion operator by calcu- 
lating the eigenvalues {—Z u } and the corresponding eigen-functions {|?/v)}. Hereafter we 
use Greek indices to denote the eigenstates and Latin indices to denote the basis states of 
the non-adiabatic diffusion operator. Because the non-adiabatic Liouville operator is non- 
Hermitian, the eigenvalues are generally given by complex values, and the right and left 
eigen-functions corresponding to the same eigenvalue are not simply the Hermitian conju- 
gate to each other. H For a given eigen- value Z u , the right and left eigen-functions of the 
non-adiabatic diffusion operator are obtained from 

C\tf) = -Z u \tf), (18) 

{ri\c = -zM\- (19) 

The method of eigenfunction solution is well known for the diffusion process on the har- 
monic potential energy surface.il For a single quadratic potential U(x) = ^muj 2 x 2 , the 
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one-dimensional Fokker-Planck operator Cfp = -^(Jnr + ft^U') can be transformed into 
the quantum mechanical Hamiltonian in imaginary time, 



with 7 = Dmuj 2 /}zbT . Since the transformed potential in Eq. fl2~T|) is just the same form as 
for a simple harmonic oscillator with zero point energy compensation, the eigenvalues and 
the eigen-f unctions for the original Fokker-Planck operator can be constructed immediately 
from the eigen-solutions of the harmonic oscillator Hamiltonian. Unlike the diffusion problem 
on the single potential energy surface, there have been limited studies on the non-adiabatic 
diffusion problem involving more than one potential energy surface. In this aspect, Cukier 
and co-workers have calculated the electron transfer rate by calculating the lowest eigenvalue 
of the non-adiabatic diffusion equation; however, their calculation was limited to the weak- 
coupling regime where the Zusman rate is applicable.il 

An important issue in solving the non-adiabatic diffusion equation for electron transfer 
is the choice of the basis functions since three different free energy surfaces are involved in 
Eq. (||): two diabatic surfaces for the population density matrix elements and one averaged 
surface for the coherence density matrix element. In this paper, the eigen-functions of £12 
are used as our basis set to represent the non-adiabatic diffusion equation. In principle, 
one could have chosen the eigen-functions of C\ or £2 as basis functions, however, in that 
case one has to evaluate appropriate Franck- Condon factors when calculating the coupling 
matrix elements even with the Condon approximation. The Fokker-Planck operator £ 12 is 
defined on the averaged harmonic potential centered at E = 0, and its eigen-solutions are 




(20) 



where \i 



2D, and the quadratic potential is 






(22) 
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(rt\£ l2 = -nQ(rt\, (23) 
where the right and left eigen functions are 

rfw= (2^4*1)* ^ (24) 

and 

= <^i^ H " hk) ■ (25) 

where H n is the nth order Hermite polynomial. As shown below, this choice of the basis set 
is convenient for our purpose. 

To be consistent with the Cu basis set, we separate the real and imaginary parts of the 

coherence density matrix, namely, u =Repi2 and v =Imp 12 , and rewrite Eq. (|9|) as 

Pl = (£ 12 + SC) Pl - 2Vv, (26a) 

p 2 = (£12 - 8C) P 2 + 2Vv } (26b) 

u = C 12 u + ui 2 v, (26c) 

v = £ 12 v - LU 12 u + V(pi - p 2 ), (26d) 



where we have defined 5C as 



6C = Al ~ £22 . (27) 



Then, all the relevant operators in Eqs. ( p6a| )-( |26o T [ ) can be evaluated in terms of the right 
and left eigen-f unctions of £ i2 , giving 

(<^|£ 12 |0£> = -n{l5 nm , (28) 



%\6£\<l£) = ^V^+lS n , m+1 , (29) 



A 



^nl^ial^m) = v2A/c jB T( v /m5 nim _i + \Jm + 15 n ,m+i) - e^nm, (30) 

(0nl^l0m) = ^nm, (31) 
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where we assume the Condon approximation, i.e., the electronic coupling matrix element 
is independent of the solvent degrees of freedom. With the basis set, we can expand the 
density matrix elements as 



Pl (E,t) = Y,a n (t)<f>«{E), 

n=0 

oo 

n=0 

oo 

n=0 

oo 

v(E,t) = J2d n (t)<f>*(E). 



(32a) 
(32b) 
(32c) 
(32d) 



n=0 



Substituting Eqs. (|32a|) - (|32d|) into the eigenvalue equation Eq. (|Tq) , we have the following 
coupled linear equations 



Z v a Ti 



-Z u b n 



-Z v c n 



-Z v d n 



-nOa v — 0\ 



-nOb n + 0\ 



A 



2k B T 



\/na n _i — 2Vd n 



A 



2k B T 



nQ.Cn + J2Xk B T(Vn + ld n+1 + y/nd n -i) - ed n , 



nOd n - J2Xk B T(Vn + lc n+1 + ^/nc n -i) + ec n + V(a n - b n ), 



(33a) 

(33b) 
(33c) 
(33d) 



which is an explicit basis set representation for the two-state diffusion operator in Eq. (^j). 
The linear equations for the left eigen-solution as defined by Eq. (|l9|) can be written by the 
transpose of Eqs. ( p3a| )-( 33d ). Diagonalizing the 4N x 4N matrix (N = number of basis 
functions) defined in Eqs. ( ^3a| )-( P3^D , we obtain the eigenvalues — Z u and the corresponding 
eigenvectors of the non-adiabatic diffusion operator, 



'/'.\0nl; 



(34) 
(35) 



where R nt , and L v „ are elements of the transformation matrices. 
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In general, due to the non-Hermitian nature of the non-adiabatic diffusion operator, the 
right and left eigen-functions do not form an orthogonal set by themselves. However, when 
the eigenvalues are all non-degenerate, the left and right eigen-functions form an orthogonal 
and complete set in dual Hilbert space.E] Explicitly, we have 

= (36) 

n=0 

for the orthogonality and 

^ ] Rnv-Lvm &nmi (37) 
v 

for the completeness. Using these properties, we can construct the real time propagator for 
the operator C as 

G(t) = £lt^><V£|e-*' t , (38) 

V 

and express the time evolution of the density matrix by projecting a given initial distribution 
onto the eigenstates, giving 

\p(t)) = G{t)\p(0)) = £ \^){^W))e- Zvt - (39) 

V 

Hence, the eigen-solution to the two-state non-adiabatic diffusion equation leads to a com- 
plete description of electron transfer dynamics. 

III. RESULTS AND DISCUSSIONS 

In the section, we present the spectral structure of the non-adiabatic diffusion operator 
by diagonalizing its matrix representation in Eqs. (|33a|)-(]3"3"ciD. In principle, we need infinite 



number of basis functions to diagonalize the non-adiabatic diffusion operator, however, in 
practice, we have to truncate our basis set at some finite number. In all the calculations 
below, we have used iV = 50 — 200 to diagonalize the 4N x AN matrix and the effect of finite 
number basis on the spectral structure has been carefully examined. 
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A. Spectral Structure 



1. Mixed-valence systems 

In the mixed-valence compounds, the electronic coupling constant has the same order 
of magnitude as the reorganization energy and the electron transfer dynamics is usually 
probed experimentally at room temperature in polar solvents. To study this process, Evans, 
Nitzan, and Ratnerl carried out real time path-integral simulations for the photo-induced 
electron transfer reaction in (NH3)5Fe n (CN)Ru m (CN)5. Based on their model, we chose the 
parameters for the calculation shown in Fig. 1 as /3Q = 0.6716, (3\ = 18.225, f3V = 11.99, 
and (3e = 18.705. As mentioned in the introduction the mapping between the spin-boson 
Hamiltonian and the semi-classical Zusman equation is not rigorously defined. For example, 
for the non-adiabatic diffusion equation, the solvation energy correlation function takes 
an exponential form with the rate Q, whereas, for the spin-boson model Hamiltonian, it 
depends on the functional form of the spectral density. It can be shown that the Ohmic 
spectral density with an exponential cut-off uj c 

J(uj) = i]uexp(— u/lj c ), (40) 

used in the calculation of Evans et ai, leads to an energy correlation function with a 
Lorentzian form at high temperature,^ 

2r]u c k B T 1 

C SB (t) « — — - 41 

7r i + {uj c ty 

Then, the relaxation rate Q used in our calculation is taken as the inverse of the mean 
survival time of CsB{t), which is Q = 2uj c /tt. 

In Fig. 1 the spectral structure of the non-adiabatic operator is shown in complex space. 
We have used A = 200 (4 A = 800) basis functions to calculate the eigenvalues. To remove 
the effect of finite basis set from the resulting spectral structure, we only show the first 

13 



400 eigenvalues in the complex plane. Since the non-adiabatic diffusion operator is non- 
Hermitian, the resulting spectrum shows complex conjugate paired eigenvalues as well as 
real eigenvalues, giving rise to the tree structure with three major branches (which we will 
call the eigen-tree) . In Fig .1, we separate the real and imaginary parts of eigenvalue by 

— Z v — -k v - iuj v . (42) 

Obviously, the real part, k u , is always negative as all non-equilibrium physical quantities 
decay to zero at time infinity, and it scales linearly with the index v since the relaxation rate 
corresponding to the nth basis state 0„ is proportional to n. In general, the relative magni- 
tudes of real and imaginary parts of eigenvalues determine the time-evolution of the density 
matrix: the real eigenvalues correspond to the simple exponential decay components and 
the complex conjugate paired eigenvalues correspond to the damped oscillation components. 

To classify the eigenvalues quantitatively according to their dynamic behavior, we intro- 
duce the dimensionless quantity 6 V 

0, = ^f, (43) 

where k v is the decay rate and 2n/uj u is the oscillation period. The time-evolution of the den- 
sity matrix component associated with the eigenvalue Z v is an exponential decay if 6 U = oo, 
an under-damped oscillation if 9 U > 1, and a damped oscillation if 9 V < 1. The relative 
amplitude of the each component depends on the overlap matrix element between the ini- 
tial density matrix and the eigenstate. As an approximate criterion for the classification 
of the eigenvalues, the slope corresponding to V — 1 is shown in the eigen-tree diagram 
in Fig. 1. There are a few eigenstates around and below the 9 V — 1 line, with a typical 
rate of (3k v as 5. For the parameters used in the calculation, (3 corresponds to ~170 fs in 
real time, and, therefore, these eigenstates exhibits damped oscillations with a period and 
a decay time in the femtosecond regime. In their real-time path integral simulations, Evans 
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et al. showed that the population in the acceptor state oscillates with a few femtosecond 
period and these oscillation decays in within 20 femtoseconds. Thus, qualitative features 
of the electron transfer dynamics can be predicted and understood from a careful exami- 
nation of the spectral structure. Since the spectral analysis presented here is based on the 
semi-classical diffusion equation while the path-integral study is based on the quantum me- 
chanical spin-boson Hamiltonian, the comparison between the two approaches is expected 
to be qualitative. In the following subsection, further analysis reveals the nature of these 
oscillations. 

2. Dependence on the coupling constant V 

To examine the underlying spectral structure in more details, eigenvalues of the non- 
adiabatic diffusion operator are plotted as functions of the electronic coupling constant in 
Fig. 2. All the parameters except for the electronic coupling constant are the same as used 
in Fig. 1. 

In Fig. 2(a), the real parts of the first 20 eigenvalues are shown as functions of the 
electronic coupling constant. Note that eigenvalues corresponding to complex conjugate 
pairs have the same real part, thus they coalesce in the real eigenvalue diagram. When 
the coupling constant is very small (f3V <C 1), the real part of the first non-zero eigenvalue 
is very well separated from the eigenvalues of excited states, so the dynamics of electron 
transfer can be considered as a incoherent rate process with a well-defined rate constant, 
k\. When the coupling constant is larger (/3V ~ 1), the first excited state becomes close 
to the second excited state, and they start to merge into a complex conjugate pair. If the 
coupling constant increases further, eigen-values show a bifurcation behavior at f3V 10. 
Therefore, in this regime, the electron transfer kinetics show multiple time-scale relaxation 
as well as coherent oscillation. The complicated behavior of coalescence and bifurcation in 
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the real eigenvalue appears more frequently at higher states. 

Another interesting feature of the real eigenvalue diagram is that a set of real eigenval- 
ues decreases consistently as the coupling constant increases from zero. It turns out that 
these eigenstates take on large imaginary parts, which are responsible for the onset of the 
imaginary branches of the eigen-tree. In Fig. 2(b), the imaginary parts of the lowest 30 
eigenvalues are plotted as functions of the coupling constant. Interestingly, the imaginary 
part of the eigenvalue increases approximately linearly with the coupling constant at large 
coupling regime. In fact, the dependence on the coupling constant is similar to that of the 
Rabi frequency for the two-level system, 

n R = Ve 2 + AV 2 , (44) 

which is shown in Fig. 2(b). As pointed out in a recent paper ,@ electronic coherence in 
mixed-valence systems arises from Rabi oscillations between two adiabatic surfaces and 
decays because of dephasing. 

To demonstrate the correlation of the real and imaginary parts of the eigenvalues as 
functions of the coupling constant, we present a three dimensional plot of the spectral 
structure in Fig. 2(c). For clarity, only the positive branches of the imaginary eigenvalues 
are shown. If we compare Fig. 2(c) with Fig. 2(a), the very rapidly decaying states shown in 
Fig. 2(a) take on large imaginary parts corresponding to the Rabi oscillations as the coupling 
constant increases, and these states are responsible for the onset of the imaginary branches 
in the eigen-tree for the mixed- valence system shown in Fig. 1. 

B. Density Matrix Propagation 

To check the validity of the spectral analysis as a density matrix propagation scheme, 
we calculated the time-evolution of the density matrix by applying the propagator defined 
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by Eq. (|38|) to the initial density matrix for various energy biases. Although it may seem 
straightforward to use the spectral method as a propagation scheme, the case for a non- 
Hermitian operator is not trivial and has not been explored. The main reason is that 
though the left and right eigen-functions of a non-Hermitian operator can be shown to form 
a bi-orthogonal set for the non-degenerate eigenvalue case, numerically these eigen-functions 
may not be stable enough to be used as a complete orthonormal basis for the density matrix 
propagation, especially in the nearly degenerate eigenvalue case. We can understand the 
situation as follows: When the two nearly degenerate eigenvalues Zi and Z 2 are obtained 
from a non-Hermitian operator, the orthogonality implies that (L 2 \ and \R\) are orthogonal 
to each other as well as (L\\ and |i?2)- When two eigenvalues become very close to each 
other, unlike the Hermitian operator case, (Li| and (L 2 \ almost coincide and so do \Ri) 
and I-R2), so that (L\\ and become almost orthogonal to each other. To still satisfy 
the normalization condition (L n \R n ) = 1 in this case, the eigenfunction should be scaled 
up, thus making the spectral structure very sensitive to the numerical error involved in the 
calculation of eigenfunctions. For an interesting discussion on this point, one may refer 
to the work by Nelson and co-workers.il Due to this numerical instability, the use of the 
spectral method as a density matrix propagation scheme is not without limitation. 

Figure 3a shows the spectral structure and the time-evolution of the density matrix 
propagation for the case of /3Q = 1, (3X = 15, f3V = 12, and j3e = 5. Generally, when 
the energy bias is small {(3e < 5), the left and right eigenfunctions can form a complete 
orthonormal basis set, so the spectral method is stable and can be used as a numerical 
propagation method for the density matrix. With a large energy bias, however, the calculated 
eigenfunctions may not form a complete orthonormal basis. To model for the photo-induced 
back electron transfer experiment in the mixed-valence compounds the initial density matrix 
is chosen as a thermal equilibrium distribution of the donor state(i.e. 1-state) pumped to 
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the acceptor state(i.e. 2-state),ilI@ 




(g + A) 2 
2A| 



) 



4 



(45a) 



p 12 (£,0) =p 21 (£,0) = 0. 



(45b) 



It would be straightforward to calculate the spatial distribution of the density matrix in time 
p(E, t) by applying the propagator in Eq. fl38|) to the above initial density matrix; however, 
to demonstrate the overall temporal behavior only the time evolution of the total population 
in the acceptor state is calculated, 



In order to check the validity of the spectral method as a propagation scheme in this case, 
we also calculated the time evolution of the density matrix by directly solving the AN 
differential equations for the expansion coefficients of the density matrix using the Bulirsh- 
Stoer algorithm,^ and the comparison in Fig. 3(a) shows a perfect agreement. If only the 
transient behavior is concerned with, the direct propagation method would be preferred over 
the spectral method, however, the spectral propagation has the advantage when calculating 
the long time behavior once the complete spectrum is known. Overall, the computational 
costs for two method are comparable to each other. As expected from the spectral structure 
shown in the previous section the population in the acceptor state shows an underdamped 
coherent oscillation behavior at initial times followed by a damped oscillation behavior at 
later times. 

Further, we have also studied the density matrix propagation for different energy biases 
to examine the electronic dephasing effect. As seen from Fig. 4(a), the increase in energy 
bias destroys the electronic coherence dramatically. Another interesting observation is the 
phase shift in the population dynamics as the energy bias is varied, and it is because the 
Rabi oscillation frequency increases with energy bias. We can confirm the temporal behavior 




(46) 
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of the density matrix propagation by examining the spectral structure shown in Fig. 4(b). 
The period of the initial coherence is estimated to be t osc m 0.25/5 from Fig. 4(a). In 
comparison, the Rabi frequency for the corresponding adiabatic two-level system is given by 
Qr ■ = Ve 2 + 4V 2 ~ 25/3" 1 , which can also be obtained from the onset of imaginary branches 
in the eigen-tree shown in Fig. 4(b), and the estimation is consistent with the oscillation 
period observed in the dynamics since r osc m 2tt/Qr. The real eigenvalues of the lowest 
excited states in the the imaginary branches are estimated to be f3k m 1 — 2, and they agree 
with the decay time of the oscillation amplitude in Fig. 4(a), confirming the validity of the 
spectral method as a density matrix propagation scheme. Even though it has been well 
known in the literature that the damping of population is enhanced with increased energy 
asymmetry,0 we have also confirmed this through the spectral analysis method. 

As an example of the eigenfunction responsible for the coherent oscillation behavior ob- 
served in Fig. 4 (b), we show the left and right eigenfunctions corresponding to a complex 
eigenvalue (3Z = 2.6228 ± z26.394 for a symmetric case (fie = 0) and (3Z = 2.8057 ± z26.466 
for an asymmetric case(/3e = 5) in Figs. 5 and 6. The eigenfunctions corresponding to 
a complex conjugate pair of eigenvalues are also complex conjugate to each other; there- 
fore, the frequency spectrum of the density matrix evolution is proportional to the norm of 
wavefunction. We note that the left eigenfunction is more extended than the right eigenfunc- 
tion. Although the population distribution in the donor and acceptor states corresponding 
to coherent oscillation is inverted with respect to the Boltzmann distribution, it does not 
contribute to the steady-state population distribution due to the transient nature. 

IV. CONCLUDING REMARKS 

In this paper we have applied the spectral analysis method to the non-adiabatic two-state 
diffusion equation, that describes electron transfer dynamics in Debye solvents. In particular, 
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we have examined electronic coherence in mixed-valence compounds, and demonstrated 
that underdamped Rabi oscillations are observed in an overdamped solvent environment. 
Detailed study of the spectral structure of the non-adiabatic operator for various energy 
biases and coupling constants allows us to determine the underlying mechanisms of electron 
transfer kinetics. Eigenvalues form three branches in the eigen-diagram: a single branch of 
real eigenvalues and two symmetric branches of complex conjugate eigenvalues. In strongly 
coupled systems, all three branches have a similar order of magnitude, indicating that both 
multiple-exponential decay and coherent oscillations can be observed experimentally. 

We have investigated the dependence of the spectral structure on the coupling con- 
stant. In the very weak coupling regime, the lowest excited state is well separated from 
higher states, which makes the electron transfer dynamics a well-defined rate process. In 
the strong coupling regime, however, the eigenvalue diagram shows coalescence/bifurcation 
behavior in the complex plane. We have used the spectral method to calculate the time- 
evolution of the density matrix, and indeed, observed electronic coherence in the temporal 
behavior of population in the acceptor state for non-equilibrium initial distributions. We 
also found a good agreement between results of the spectral propagation method and of 
the numerical propagation method for small energy bias cases. Due to non- Hermit ianity of 
the non-adiabatic operator, the spectral propagation method was not numerically stable for 
large energy bias cases. 

For an isolated quantum system, the eigen-solution to the Schrodinger equation com- 
pletely determines its dynamics. In a similar fashion, the eigen-solution to the non-adiabatic 
diffusion operator completely characterizes the dynamics of a dissipative system and thus 
provides a powerful tool to analyze dissipative dynamics. It is well known that quantum dy- 
namics comes from the underlying spectra, especially in gas-phase chemical systems;!! how- 
ever, the spectral aspect of condensed phase dissipative systems has not been well recognized 
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yet and deserves further investigation. Though the analysis presented here is restricted to 
semi-classical dissipative systems, it may also be applied to quantum dissipative dynamics. 
In principle, we can derive the evolution equation for quantum dissipative systems either 
from first principles or through numerical reduction, and then pose the quantum dissipative 
equation of motion as an eigen-value problem. Along this line, the dissipative dynamics of 
the spin-boson Hamiltonian, which has been studied mostly as a dynamic problem,!'^ can 
also be explored spectral problem in the future. 
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FIGURES 

FIG. 1. A plot of the lowest 400 eigenvalues for the non-adiabatic operator in a mixed- valence 
system. The parameters are : (59, = 0.6716, j3\ = 18.225, j3V = 11.99, and /3e = 18.705. The 
dot-dashed line is for the case k = lo/2tt. 

FIG. 2. Plots of (a) real and (b) imaginary parts of the lowest 30 eigenvalues as a function of 
the coupling constant, V. Except for the coupling constant, all the other parameters are set equal to 
those used in Fig. 1. In Fig. 2(b), open circles correspond to the Rabi frequency VLr = \/e 2 + 4V 2 . 
Figure 3(c) shows the three dimensional plot of eigenvalues as a function of the coupling constant. 

FIG. 3. Comparison between the result of direct numerical propagation and spectral propa- 
gation. The parameters are chosen as (30. = 1, f3\ = 15, (3V = 12, and fie = 3. 

FIG. 4. Comparison of (a) the dynamics and (b) the spectra in the mixed-valence system for 
three different energy biases. Except for the energy bias, all the other parameters are set equal to 
those used in Fig. 3. Agreements between the results of numerical and spectral propagation have 
been checked in these cases. 

FIG. 5. (a) Right and (b) left eigenfunctions with an eigenvalue [3Z = 2.6228 ± i26.394 for a 
symmetric bias case.(/3e = 0) All the other parameters are set equal to those used in Fig. 3 except 
for the energy bias. Each line corresponds to p\ (solid), p2 (dashed), u(dot-dashed), and v (dotted), 
respectively. 

FIG. 6. (a) Right and (b) left eigenfunctions with an eigenvalue [3Z = 2.8057 ± i26.466 for an 
asymmetric bias case.(/3e = 5) All the other parameters are set equal to those used in Fig. 3 except 
for the energy bias. Each line corresponds to p\ (solid), p2 (dashed), u(dot-dashed), and v (dotted), 
respectively. 
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